\documentclass[11pt]{beamer}
\usepackage[UTF8,scheme=plain]{ctex}
\usepackage{listings}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{graphicx}
\usetheme{Boadilla}

\usepackage{framed} % 可以用 \begin{shaded}，即背景色块
\definecolor{shadecolor}{rgb}{0.9,0.9,0.9}

\lstset{
	basicstyle=\footnotesize,
    flexiblecolumns,
    showstringspaces = false,
	numbers=left,
    numberstyle=\tiny,
    numbersep=4pt,
	breaklines,%自动换行
	frame=shadowbox,%加边框
	backgroundcolor = \color{red!30!green!30!blue!30},
	%背景色：蓝灰色
	language=C++, % 语言为 c++
	commentstyle=\color{red!80!green!80!blue!80},% 注释用蓝灰色字体标出
	keywordstyle=\color{blue!70}, % 高亮部分用蓝色字体
	tabsize = 4,
    escapeinside=@@
}

\begin{document}
	\author{ 路毅 \hspace{0.3cm} 曲阜师范大学 }
	\date{\number\year 年 \number\month 月 \number\day 日}
	\title{计算物理：一维波动方程求解}
	%\subtitle{}
	%\logo{}
	%\institute{}
	%\date{}
	%\subject{}
	%\setbeamercovered{transparent}
	%\setbeamertemplate{navigation symbols}{}
	\begin{frame}[plain]
	\maketitle
\end{frame}

\kaishu

\begin{frame}[fragile]{1维波动方程}
一维波动方程（比如绳波）的方程为
\begin{eqnarray}
\frac{ \partial^2 }{ \partial x^2 } u  = \frac{ \partial^2 }{ \partial t^2 }u,
\nonumber
\end{eqnarray}
式中$u(x,t)$代表绳子上$x$处质点在时刻$t$的位移。
\vspace*{\fill}
\end{frame}

\begin{frame}[fragile]{1维波动方程}
一维波动方程（比如绳波）的方程为
\begin{eqnarray}
\frac{ \partial^2 }{ \partial x^2 } u  = \frac{ \partial^2 }{ \partial t^2 }u,
\nonumber
\end{eqnarray}
式中$u(x,t)$代表绳子上$x$处质点在时刻$t$的位移。

如果初始时1m长的绳子没有形变，各处静止，则初始条件为
\begin{eqnarray}
u(x,t=0) &=& 0, \nonumber \\
\frac{\partial}{\partial t} u |_{t=0} = 0. \nonumber
\end{eqnarray}
\vspace*{\fill}
\end{frame}

\begin{frame}[fragile]{1维波动方程}
一维波动方程（比如绳波）的方程为
\begin{eqnarray}
\frac{ \partial^2 }{ \partial x^2 } u  = \frac{ \partial^2 }{ \partial t^2 }u,
\nonumber
\end{eqnarray}
如果初始时1m长的绳子没有形变，各处静止，则初始条件为
\begin{eqnarray}
u(x,t=0) &=& 0, \nonumber \\
\frac{\partial}{\partial t} u |_{t=0} = 0.
\nonumber
\end{eqnarray}
绳子一端做简谐振动，另一端固定在墙上，边界条件为
\begin{eqnarray}
u(0,t) &=& sin(t), \nonumber \\
u(1,t) &=& 0. \nonumber
\end{eqnarray}
{\bf 有了这些条件，波动方程的解是唯一的。}
\end{frame}

\begin{frame}[fragile]{差分形式：２阶导数}
函数的2阶导数可以用差分公式来近似计算，
\begin{eqnarray}
f''(a) \approx \frac{ f(a+h) + f(a-h) - 2 f(a) }{h^2} + O(h^2), \nonumber
\end{eqnarray}

\end{frame}

\begin{frame}[fragile]{差分形式：２阶导数}
函数的2阶导数可以用差分公式来近似计算，
\begin{eqnarray}
f''(a) \approx \frac{ f(a+h) + f(a-h) - 2 f(a) }{h^2} + O(h^2), \nonumber
\end{eqnarray}
如果把微分方程用有限差分代替，即
\begin{eqnarray}
\frac{ \partial^2 }{ \partial x^2 } u  = \frac{ \partial^2 }{ \partial t^2 }u,
\nonumber
\end{eqnarray}
变为：
\begin{eqnarray}
\frac{ u_{i+1,j} + u_{i-1,j} - 2 u_{i,j} }{(\Delta x)^2} + O((\Delta x)^2)
= \frac{ u_{i,j+1} + u_{i,j-1} - 2 u_{i,j} }{(\Delta t)^2} + O((\Delta t)^2),
\nonumber
\end{eqnarray}
其中$u_{i,j} \approx u(i*\Delta x , j*\Delta t)$。
\end{frame}

\begin{frame}[fragile]{数值递推式}
如果取$\alpha = \frac{(\Delta t)^2}{(\Delta x)^2}$，就得到显式的递推式
\begin{equation}
{\color{blue}u_{i,j+1} = (2-2\alpha) u_{i,j} + \alpha u_{i+1,j} + \alpha u_{i-1,j} - u_{i,j-1}. }
\nonumber
\end{equation}
只要知道$u_{*,j-1}$与$u_{*,j}$，就可以递推得出$u_{*,j+1}$。
\end{frame}

\begin{frame}[fragile]{数值递推式}
如果取$\alpha = \frac{(\Delta t)^2}{(\Delta x)^2}$，就得到显式的递推式
\begin{equation}
{\color{blue}u_{i,j+1} = (2-2\alpha) u_{i,j} + \alpha u_{i+1,j} + \alpha u_{i-1,j} - u_{i,j-1}. }
\nonumber
\end{equation}
只要知道$u_{*,j-1}$与$u_{*,j}$，就可以递推得出$u_{*,j+1}$。

第一步前推时，我们只知道$u_{*,0}$，不知道$u_{*,-1}$，那么可以利用初始条件
\begin{equation}
\frac{\partial}{\partial t} u |_{t=0} = 0 \Rightarrow 
u_{i,-1} = u_{i,1}, \nonumber
\end{equation}
带入递推式，得到
\begin{equation}
{\color{blue}u_{i,1} = (1-\alpha) u_{i,0} + \alpha \frac{u_{i+1,0} + u_{i-1,0}}{2}. } \nonumber
\end{equation}
这样就可以完成前推了。
\end{frame}

\begin{frame}[fragile]{代码实现}
下面的代码中是一个函数，在给定初始条件，边界条件下，计算一维波动方程，并将各个时间节点的所有$u_{*,j}$作为一行输出到文件中。
\begin{lstlisting}
void wave1d(    double xa, int nx, double dx,
                double (*boundary1)(double t), double (*boundary2)(double t),
                int nt, double dt, double (* u_initial)(double x),
                string file ){

        int i,j;
        double alpha = dt*dt/dx/dx;

        double *u_old = new double [nx]; // u_{i-1}
        double *u = new double [nx]; // u_i
        double *u_new = new double [nx]; // u_{i+1}

        u_old[0] = boundary1(0);
        u_old[nx-1] = boundary2(0);
        for(i=1;i<nx-1;i++) u_old[i] = u_initial( xa + i*dx);
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{代码实现}
\lstset{firstnumber=last}
\begin{lstlisting}
        u[0] = boundary1(dt);
        u[nx-1] = boundary2(dt);
        for(i=1;i<nx-1;i++) u[i] = (1-alpha)*u_old[i] + alpha/2*( u_old[i+1] + u_old[i-1] );

        ofstream fp(file);
        if(!fp){
                cout<<"error: failed to open "<<file<<endl;
                exit(1);
        }
        for(i=0;i<nx;i++)fp<<u_old[i]<<"\t"; fp<<endl;
        for(i=0;i<nx;i++)fp<<u[i]<<"\t"; fp<<endl;
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{代码实现}
\lstset{firstnumber=last}
\begin{lstlisting}
        for(j=2;j<=nt;j++){
                u_new[0] = boundary1( j*dt);
                u_new[nx-1] = boundary2( j*dt);
                for(i=1;i<nx-1;i++)
                        u_new[i] = (2-2*alpha)*u[i] + alpha * u[i+1] + alpha* u[i-1] - u_old[i];
                for(i=0;i<nx;i++)fp<<u_new[i]<<"\t"; fp<<endl;

                for(i=0;i<nx;i++){
                        u_old[i] = u[i];
                        u[i] = u_new[i];
                }
        }
        fp.close();
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{边界条件}
下面的代码设定初始条件和边界条件。
\begin{lstlisting}
double u_initial(double x){
        return 0;
}

double boundary1(double t){
        return sin(10*t);
}

double boundary2(double t){
        return 0;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{主函数}
\lstset{firstnumber=last}
\begin{lstlisting}
int main(){

        int i,j,nx=300;
        double xa=0, xb=3;
        double dx = (xb-xa)/nx;
        int nt=500;
        double ta=0, tb=5;
        double dt = (tb-ta)/nt;
        double **matrix;

        string file = "wave1d.txt";

        wave1d( 0, nx, dx, boundary1, boundary2,
                nt, dt, u_initial, file);
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{主函数}
\lstset{firstnumber=last}
\begin{lstlisting}
        // transpose the data matrix
        matrix = new double *[nt+1];
        for(i=0;i<nt+1;i++) matrix[i] = new double [nx];

        ifstream fp1( file );
        for(i=0;i<nt+1;i++)
                for(j=0;j<nx;j++)
                        fp1>> matrix[i][j];
        fp1.close();

        ofstream fp2( file );
        for(j=0;j<nx;j++){
                fp2<< xa + j*dx <<"\t";
                for(i=0;i<nt+1;i++)
                        fp2<< matrix[i][j]<<"\t";
                fp2<<endl;
        }
        fp2.close();

        return 0;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{作图代码}
代码中最后一段是将数据文件的行和列做转置，方便用gnuplot画成不同的桢，然后生成动图gif文件。
gnuplot代码如下
\begin{lstlisting}
set term gif animate delay 1
set output "wave1d.gif"
set xrange [0:3]
set yrange [-2:2.5]
set xlabel "x(m)"
set ylabel "y(m)"
do for[t=2:502:1]{
        plot "wave1d.txt" u 1:t w l lw 3 t sprintf("string wave: v=1m/s, t=%f s", (t-2)*0.01)
}
set output
\end{lstlisting}
生成的动图可用作演示。
\end{frame}

\begin{frame}[fragile]{总结：1.波的传播}
\begin{figure}
\includegraphics[width=0.6\textwidth]{src/wave1d_propagation}
\end{figure}
客观世界与数值世界遵循{\color{blue} 同样的数学方程}，计算机就会再现客观世界发生的{\color{blue} 同样的物理现象。}
\end{frame}

\begin{frame}[fragile]{总结：2.波的干涉、驻波}
\begin{figure}
\includegraphics[width=0.6\textwidth]{src/wave1d_standingwave}
\end{figure}
行波在右端点反弹，反射波向左传播，与向右的入射波发生干涉，形成驻波。驻波的波节静止不动，波腹的振动幅度为入射波的2倍。
\end{frame}

\begin{frame}[fragile]{总结：3.反射半波损失}
\begin{figure}
\includegraphics[width=0.6\textwidth]{src/wave1d_halfwaveloss}
\end{figure}
行波在静止的右端点反弹时，相位增加$\pi$，这称作半波损失，在光波中也存在。
另有一个动图对此进行说明。
\end{frame}

\end{document}